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The hydrodynamic interactions of a suspension of self-propelled particles are studied using a 
direct numerical simulation method which simultaneously solves for the host fluid and the swimming 
particles. A modified version of the "Smoothed Profile" method (SPM) is developed to simulate 
microswimmers as squirmers, which are spherical particles with a specified surface-tangential slip 
velocity between the particles and the fluid. This simplifled swimming model allows one to represent 
different types of propulsion (pullers and pushers) and is thus ideal to study the hydrodynamic 
interactions among swimmers. We use the SPM to study the diffusive behavior which arises due to 
the swimming motion of the particles, and show that there are two basic mechanisms responsible for 
this phenomena: the hydrodynamic interactions caused by the squirming motion of the particles, and 
the particle-particle collisions. This dual nature gives rise to two distinct time- and length- scales, 
and thus to two diffusion coefficients, which we obtain by a suitable analysis of the swimming motion. 
We show that the collisions between swimmers can be interpreted in terms of binary collisions, in 
which the effective collision radius is reduced due to the collision dynamics of swimming particles in 
viscous fiuids. At short time-scales, the dynamics of the swimmer is analogous to that of an inert 
tracer particle in a swimming suspension, in which the diffusive motion is caused by fluid-particle 
collisions. Our results, along with the simulation method we have introduced, will allow us to gain 
a better understanding of the complex hydrodynamic interactions of self-propelled swimmers. 
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I. INTRODUCTION 



Swimming microorganisms, from bacteria, to algal cells, to spermatozoa, are ubiquitous in biological processes, and 
even though the specific propulsion mechanism can vary, the motion of these microswimmers is defined by two basic 
characteristics: (1) they are moving at very low Reynolds number, where viscous forces are dominant (inertial forces 
can be neglected), and thus (2) the net force on the body (including the cilia or flagella used to generate the motion) 
is zero[Tl Although we have a clear understanding of how these organisms can generate motion, the physical 
properties of a suspension of such swimmers are still not completely understood [SHZ]- In particular, the rheological 
properties show a non-trivial dependence on the concentration of swimmers [SVllOj . which deviates considerably from 
that of inert force-free colloidal particles [llj. In addition, hydrodynamic interactions among swimmers have also been 
shown to give rise to collective motion [T^l [13] • However, theoretical studies on such phenomena have either neglected 
the role of hydrodynamic interactions, or used only far-field approximations [71 114j. 

In addition to the biological interest of understanding the hydrodynamic interactions of swimming microorganisms, 
recent experimental work has shown how to design micro-motors which can swim in the absence of external fields 
thanks to interfacial phoretic effect s fTBl - fT^ . which allow the particles to transform the local chemical or thermal 
energy into mechanical energy. Furthermore, the possibility of using these systems to perform work, whether it be to 
transport cargo[5ni[lI] or separate particle mixtureSj22j, has led to many recent experimental and simulation studies. 
Although a simple theoretical framework for these swimming systems, which can be used to study the effect of the 
shape and surface properties on the swimming motion, has already been proposed[32 |^], the collective dynamics, 
and in particular the role of the hydrodynamic interactions remains an open question. In this paper, we present a 
direct numerical simulation method (DNS) for self-propelled particles which attempts to overcome this difficulty, by 
solving the equations of motion for the host fluid as well as the swimming particles. 

A detailed description of the swimming mechanism of real microorganisms remains computationally prohibitive, 
particularly if one is interested in the collective dynamics of such systems; therefore, we must look for a simple 
generic swimmer, which can be easily modeled, but which manages to reproduce the basic flow properties of actual 
swimmers, such that the hydrodynamic interactions are accurately reproduced. The "squirmer" model we use was 
introduced over 50 years ago[551 HB] to study the propulsion of spherical ciliate particles, which generate motion 
through the synchronized beating of an envelope of cilia at their surface. Fortunately, on a mesoscopic scale, the 
effect of this surface motion can be replaced by a specified slip velocity between the particle and the fluid. In 
addition, the parameters of the model can be tuned to mimic different propulsion mechanism, allowing us to study 
the hydrodynamic interactions of different (idealized) microswimmers within a unified framework. We present a 
modified "Smoothed Profile" met hod [27], previously developed to study colloidal dispersions, which is capable of 
incorporating this squirming motion by imposing the appropriate slip boundary conditions at the surface of the 
particles. The advantage of this method is that it allows us to accurately and efficiently include the many-body 
hydrodynamic interactions among the particles, and it can be easily extended to complex host fluids. 



In this work we study the effect of hydrodynamic interactions on the dynamics of a suspension of self-propelled 
particles. The swimming motion of the particles is known to give rise to a diffusive behavior [28], even in the absence 
of thermal fluctuations, and a tracer particle placed in such a system will also undergo diffusion 29, 30 . However, the 
nature of the diffusion is different in both cases, as is evidenced by the different scaling behavior with the concentration 
of swimmers [31]. In the former case, the diffusion is caused primarily by particle-particle collisions, while the latter is 
due to the hydrodynamic interactions caused by the squirming motion of the particles. By a suitable analysis of the 
particle displacements, and the decay in the velocity fluctuations, we show that both phenomena are present in the 
motion of the squirmers themselves, as should be expected. We are thus able to extract the two underlying time- and 
length-scales of the system, corresponding to the two distinct diffusive motions, just from the particle trajectories. The 
importance of this cannot be understated, as these two mechanisms (hydrodynamics and particle collisions) are known 
to be fundamental in determining the physical properties of swimming suspensions [32] : additionally, measuring tracer 
diffusion in dense suspensions experimentally can be quite challenging, and the available simulation methods for this 
are somewhat involved |31]. Finally, using basic results from kinetic theory, we analyze the collisions of the swimmers 
to derive effective collision radii, which are shown to be independent of concentration, and which depend only on the 
squirming mode of the swimmers (i.e., the strength of the pusher/puller character). Somewhat surprising is the fact 
that these effective particle sizes are considerably smaller than the actual size of the particles, which is a consequence 
of the collision dynamics of swimmers in viscous fluids. This paper is organized as follows: Section |Tl| introduces the 
mathematical model of the swimmers and the simulation method. Section |III| validates the computational method 



against known analytical results, and Section IV presents the results of our study. 
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II. THE MODEL AND METHODS 



A. Swimmer Model: Blake's Squirmers 




FIG. 1. (color online) Schematic representation of the propulsion mechanism and flow profiles of a pusher and a puller, (a) and 
(b) respectively. These swimmers can be represented using Blake's squirming model, in which the detailed propulsion mechanism 
is replaced by a specified sUp velocity at the surface of the particles, (c) and (d), for pushers and pullers, respectively. 

We consider a simple model of self-propelled spherical swimmers, originally introduced by Lighthill[25 and later 
extended by Blake [26 , which move due to a self-generated surface-tangential velocity u^. This specific mechanism 
was proposed as a model for an ideal ciliate particle, in which the synchronized beating of the cilia at the surface gives 
rise to a net motion, in the absence of any external fields. If one assumes that the displacements of this envelope of 
cilia are purely tangential, then the effective (time-averaged) slip velocity for these squirmers is given by [26] 
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where e is the squirmer's fixed swimming axis (i.e., we consider that each squirmer carries with it a fixed coordinate 
system which determines its preferred swimming direction at each instant), r is a unit vector from the particle center 
to a point on the surface, P,' is the derivative of the n-th order Lcgendrc polynomial, and Bn is the amplitude of 
the corresponding mode. Neglecting all squirming modes higher than three, i?„ — Q{n > 3), the following simple 
expression for the surface tangential velocity, as a function of the polar angle 6 ~ cos^^ {f ■ e), is obtained 
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where a = B2/B1 determines whether the swimmer is a pusher (a < 0) or a puller (a > 0). An example of the former 
are spermatozoa and most bacteria, of the latter the unicellular alga Chlamydomonas. A schematic representation 
of the flow profile generated by these two types of swimmers is given in Figure [T] Although the squirmer model we 
adopt forgoes describing the detailed propulsion mechanism, it is capable of distinguishing between pushers/pullers 
and provides an adequate approximation to the far-field flow profile generated by these swimmers. For Newtonian 
fluids, which is the only case considered here, the swimming speed U of the squirmer is determined uniquely by the 
first mode irrespective of the size of the particle, asU — 2/3i?i, while the second mode gives the strength of the 
stresslet [531 [33] . The velocity field generated by a single such squirmer, in the Stokes regime, was solved analytically 
by Ishikawa et al.f33' and is given, in the laboratory frame (fluid at rest far away from the particle), by 
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where a is the radius of the particle. Notice that for neutral swimmers (a = 0) the velocity field decays as r~^, 
while for pushers/pullers (a ^ 0) the velocity field decays as r~^. In contrast, the velocity field for a sedimenting 
particle (or a particle experiencing a net body force) decays as r~^[3S]. This will have important consequences when 
considering the hydrodynamic interactions of suspensions of swimmers. 



B. Simulation Method : Basic Equations 



We propose a direct numerical simulation (DNS) procedure to study a system of self-propelled squirmers based on 
the "Smoothed Profile" (SP) method [271|3S|, which allows one to efficiently solve both the Navier-Stokes equation 
(NS), for the fluid motion, and the Newton- Euler equations, for the colloids. This method has been successfully used 
to study the diffusion, sedimentation, electro-hydrodynamics, and rheology of colloidal dispersions in incompressible 
fluids [371 - 140] . and recent work has shown how it can be extended to treat compressible fluids within a fluctuating- 
hydrodynamics approach|41|. A detailed error analysis of the method can be found in ref. The basic idea is to 
replace the sharp boundary at the colloid/fluid interface with a diffuse interface of finite thickness (. This allows 
one to discretize the system using a fixed Cartesian coordinate grid, since the interface will always be supported by 
multiple-grid points. Although a loss of accuracy at the surface of the particle is inevitable, we can easily impose 
periodic boundary conditions (PBC), and use a Fourier spectral method to solve for the fluid equations of motion. 
The particles in the SP method are not treated as boundary conditions for the host fluid, but rather as a body 
force in the NS equation. Thus, we avoid the mesh-reconstruction problems that plague most computational fluid 
dynamics methods for systems with moving boundaries. We are aware of two alternative simulation methods that 
aim to describe these squirmer suspensions at the same level of description, the first was developed by Ramachandran 
et al.[33] using a Lattice Boltzmann (LB) model, and the second was originally introduced by Downton and Stark^i] 
within a multi-particle collision dynamics (MPC) framework, and later extended by Gotze and Gompper|45j to recover 
the correct rotational dynamics. Although the implementation details are specific to each of the models (LB, MPC, 
SP), the general mechanism used to obtain the squirming motion is the same in all three cases: local conservation of 
momentum. For the moment though, these DNS approaches have not been extensively used to study these types of 
swimming systems; the most popular methods, which still account for the hydrodynamic interactions, have usually 
been based on Stokesian D vnamics [51 [5T1 [551 H51 - H5] . and are thus limited to Newtonian fluids in the Stokes regime. 

In what follows, we briefly review the governing equations for a dispersion of inert colloids in a simple Newtonian- 
fluid, before considering how the equations must be modified for use with the SP method for swimming particles 
with slip boundary conditions. The formulation we present closely follows that of ref. I36[ in which a more detailed 
description of the computational algorithm can be found. The motion of the host fluid is determined by the Navier- 
Stokes equation with the incompressibility condition 

V • M/ = (4) 
p{dt+Uf ■V)uf = V ■ a- (5) 

where p is the total mass density of the fluid, «/ is the host fluid velocity field, and cr is the stress tensor 

(T = -pi + cr' (6) 

er' = T] Vu/ + (Vtt/)* (7) 

with rj the shear viscosity of the fluid. Consider a mono-disperse system of A^-spherical particles, of radius a, mass 
Mp, and moment of inertia Ip — 2/5Mpa'^I (with I the unit tensor). The evolution of the colloids is given by the 
Newton-Euler equations [49]. 

R^ = V, Q, = Q, skew {Vli) (8) 

where Ri and Vi denote the center of mass position and velocity of particle i, respectively, Qt is the orientation 
matrix [50] and fli the angular velocity, with skew(f2.i) the skew-symmetric angular velocity matrix 



skew(rj,) = \ ni (9) 




The forces on the particles are comprised of hydrodynamic contributions arising from the fluid-particle interactions 
F^, the colloid-colloid interactions due to the core potential of the particles F'~^ (which prevents particle overlap), and 
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a possible external field contribution F^^^ (such as gravity). Likewise, the torques on the particles can be divided into 
a hydrodynamic and an external contribution TV'^'^' (for simplicity, the particle-particle interactions are assumed 
to be given by a radial potential). In what follows we consider buoyancy- neutral particles, so that F'^^^ = TV™' — 0. 
Finally, conservation of momentum between the fluid and the particles implies the following hydrodynamic force and 
torque on the i-th particle 

F^ = JdS, (7 (10) 
ATf = J [x - R,) X (dS , ■ a) (11) 
where J dSi indicates an integral over the particle surface. 



C. Simulation Method: Smoothed Profile Squirmers 



We now present the computational algorithm used to simulate the motion of spherical particles, with a given 
surface tangential slip velocity , using the SP method. We require that all field variables be defined over the entire 
computational domain (fluid -I- particle). The concentration field for the colloids is given as 4'{x,t) = X^iLi 4>i{^T't)i 
where 0i £ [0, 1] is the smooth profile field of particle i. This field is defined such that it is unity within the particle 
domain, zero in the fluid domain, and smoothly interpolates between both within the interface regions. Details on 
the specific definition and the properties of this profile function can be found in ref. 1271 The particle velocity field is 
defined in a similar fashion, as 



N 



(12) 



with Ti = X — Ri, which allows one to define the total fiuid velocity field as 

u{x,t) = (1 — '/')«/ + (pUp 



(13) 



where the incompressibility condition is satisfied on the entire domain V • tt — 0. The evolution equation for u is then 
derived assuming momentum-conservation between fluid and particles [27] 



p{dt + u-V)u^V -0- + p(l)fp + pf, 



sq 



(14) 



where (/>/p represents the force density field needed to maintain the rigidity constraint on the particle velocity field 
and fsq is the force density field generated by the particles' squirming motion. 

We use a fractional step approach to update the total velocity field. Let u" be the field at time t„ = nh [h is the 
time interval). We first solve for the advection and hydrodynamic viscous stress terms, and propagate the particle 
positions (orientations) using the current particle velocities. This yields 



u = u 



dsV • 



{—p* I + cr') — uu 



dsV, 



t„ + h 



ds Qiskew (fti) 



(15) 
(16) 
(17) 



where the pressure term p* in Eq. ( |15[ ) is determined by the incompressibility condition V ■ u* — 0. The remaining 
updating procedure applies to the slip condition at the particle boundary and the rigidity constraint on the velocity 
field. 

We now consider the momentum change needed to maintain the slip velocity at the surface of each of the squirmers, 
where the slip profile is imposed with respect to the particle velocities {'V^';f2-}, using the previously updated 
positions and orientations {iZ""''^; Q""*"^}. We note that at this point we do not yet know the updated particle 
velocities {V"^^; fl"'^^}, which are the values that should be used when enforcing the surface slip profile V/ = V"^^ 
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{fl'f = fl'^'^^). Therefore, we adopt an iterative solution, and as an initial guess, we use the particle velocities at the 
previous time step, i.e., V/ = V^" {fl'^ = H"). The updated total velocity field is now 









+h 


U** — U* + 




ds/sq 
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(18) 
(19) 



i=l 




FIG. 2. (color online) Schematic representation of the updating scheme used to enforce the slip boundary condition at the 
surface of the squirmers. Each particle is considered to exert a force on the fluid at the surface, in order to maintain the 
specified flow profile u" (red arrows) for the squirming motion. To ensure local momentum conservation, a counter-flow is 
added within the particle domain (blue arrows). 



The second term on the right hand side of Eq. ( 19 1 imposes the slip velocity profile at the surface of each of the 
squirmers; where (pi (x (1 — 4>i) |V</)i| is a smooth surface profile function which is non-zero only within the interface 
domain of the squirmer (normalized to have a maximum value of one), and zero everywhere else. The third term 
adds a counter-flow entirely within the particle domain, in such a way that local momentum conservation is obtained. 
Assuming rigid-body motion, with velocities 5Yi and (5$!;, this requires 



j Ax (t)i {SVi + 6Vl^ X ri) = 
dx ri X (f>i {5Vi + Sfli X 



= - dxLp^ {V/ + X n + ul - u*) 



dx ri X ipi (V/ + fl'^x ri + uf — u*) 



(20) 
(21) 



from which we can easily obtain the counter-flow terms SVi (Si^^i) from the particle velocities (JIQ. A schematic 
representation of this procedure, used to enforce the specific slip-boundary conditions for our model squirmers, is 
shown in Figure [2j Finally, the pressure term due to the squirming motion p^q is obtained from the incompressibility 
condition V • u** — 0. At this point, the momentum conservation is solved for the total velocity field. 

The hydrodynamic force and torque exerted by the fluid on the colloids (which includes all contributions due to the 
squirming motion) is again derived by assuming momentum conservation. The time integrated hydrodynamic force 
and torque over a period h are equal to the momentum exchange over the particle domain 



I d,s(ArH + Ar^q) 



= dx t 



J dx X 



(22) 
(23) 
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From this and any other forces on the colloids, the particles velocities are updated as 



yn+l ^ yn _^ 



t„ + h 



ds (AT" + iVf^' 



m: 



t„+h 



ds {Fp + Fr') 



t„+h 



(24) 
(25) 



We recall that we have imposed the slip profile u'^ with respect to the primed velocities {V/;rt'^}, which need not 
be equal to the final velocities of the particle at step n+ 1. To maintain consistency, we iterate over Eqs. (18l-(25) 



until convergence in the velocities is achieved. Finally, the resulting particle velocity field (/)""'"^Mp+^ is enforced on 
the total velocity field as 



= u** 



t„+h 



dscjjfp 



t„ + h 



ds (j)fp 



P 



(26) 
(27) 



with the pressure due to the rigidity constraint obtained from the incompressibility condition V • — 0. The total 
pressure field is then given hy p ^ p* + Pp + p^q- 

D. Choosing the Relevant Reference Frame 

Although we have not included thermal fluctuations in our system, the swimming motion of the particles is known 
to give rise to diffusive behavior at sufRciently high particle concentrations |25j. However, due to the self-propelled 
nature of the particles, the largest contribution to their displacement will naturally come from the their inherent 
swimming. As such, one must wait a very long time before the particles exhibit any type of diffusive behavior, and 
even then, it is difficult to establish what role the hydrodynamic interactions among neighboring particles are playing. 
All particles will be swimming in the flow field generated by their neighbors, and the interactions among them can 
give rise to motion perpendicular to the particle's preferred swimming direction, as well as provide a momentary 
impulse that can increase/decrease the velocity parallel to the swimming axis, or even change its orientation in space. 
In order to better understand this phenomena, we analyze the particle motion with respect to the frame of reference 
of the moving squirmers, as was proposed by Han et al.|51j to study the Brownian motion of ellipsoidal particles. 
We begin by decomposing the particle trajectories in terms of displacements between successive time intervals, which 
we take to be equally spaced. If R{tn) denotes the position of a tagged particle at time step n (i„ = n5t, with 5t a 
suitably small time interval), we can express the time evolution of its position {i2(ii)} as 



{R{tQ),R{to) + Ai, . . . , Rita) + ^ A,, . . . , i2(to) + ^ A„} 
where A^ — R(ti) — i?(ri_i). Using this notation, the translational diffusivity is given as |52j 

Dritn) = -^^{iRitn) - R{to)f 



(28) 



(29) 
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To study the coupling between translation and rotation of the particles, we will also consider the rotational diffusivity, 
defined as 



i9(i„)= f"dsn(s) 



(30) 
(31) 



to 



FIG. 3. Analysis of the particle displacements with respect to the particle's moving coordinate frame. The particle displacement 
A during a given time interval is decomposed into its components parallel A" and perpendicular A^ to the swimming axis. For 
a single isolated squirmer A^ — for any given time interval; for a suspension of swimmers the flow induced by the neighboring 
particles gives rise to non-zero perpendicular displacements. 



where i9 is the unbounded rotational displacement of the particle. 

In order to separate the "random" contribution given by the surrounding configuration from the particle's own 
swimming motion, we consider the displacements within the body frame of reference, 



(32) 



where tildes are used to denote quantities with respect to the coordinate-system attached to each of the squirmers 
(61,62,63), with the 63 the preferential swimming axis (see Figurejs]). The advantage of this approach is illustrated 
in Figure [4] where the trajectory of a single particle, from a suspension of pushers a = +2 at (/) = 0.1, is given in both 
representations: with respect to the lab- and body-space displacements. One can immediately see where the difficulty 
in analyzing the particle motion comes from, as the length scales for the directed (swimming) and fluctuating motion 
differ by an order of magnitude. We define an additional effective hydrodynamic diffusivity Z?^, in terms of these 
perpendicular and parallel displacements, as 




nStU 



DUtn) = 

Dritn) = ^ (25|,(t„) + 5^(i«)) 



(33) 

(34) 
(35) 



where the second term inside square brackets in Eq. ( 34 1 is required to remove the (average) contribution to the 
particle displacements from the inherent swimming motion ([/ = (V ■ 63) is the average velocity along the particle's 
swimming axis). These diffusivities provide a direct measure of the strength of the hydrodynamic interactions among 
the squirmers. The corresponding diffusion coefficients can be obtained from the long-time limit (assuming a plateau 
has been reached and the limit exists) 



Dt = hm Drit) 
Dr = lim DR{t) 

Dt = hm Drit) 



(36) 
(37) 
(38) 
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Since the parallel and perpendicular effective diffusion coefficients, and Dj^, exhibit no deviation from isotropic 
behavior (although this could be expected to change if persistent long-range order appears), we only consider the 
average effective diffusion coefficient Dt- A similar analysis can be performed for the rotational diffusivities, in order 
to study the rotations about the three particle axes independently, but it provides no additional information and will 
not be presented here. 




FIG. 4. (color online) Trajectory of a single particle in a suspension of pullers a = +2 at (jj = 0.1. (a) The real lab-space 
trajectory of the particle, and (b) the body-space trajectory. The latter is constructed by transforming the individual particle 
displacements to the frame of reference of the particle, and expressing them in terms of displacements parallel and perpendicular 
to the particle's swimming axis eg. In addition, the (parallel) motion due to the average swimming speed, Ax{t) = {V ■ e.^,)t, 
has been removed. Thus, Figure (b) gives the random motion induced on the particle by the local fluctuations in the surrounding 
flow field. The length scale is given by the particle radius a/A = 5 and the projections of the trajectories onto the bottom 
plane have been added as a visual guide. 

Finally, we also consider the velocity auto-correlation functions [52] 

Cv{t)^{V{t)-V{to)) (39) 
Cn{t) = {n{t) ■ fl(to)) (40) 

and in particular, the correlation functions for the velocity components parallel and perpendicular to the swimming 



C\,{t) = {V\\{t)-VHto)) (41) 
cl{t)^(vHt)-V\\{t^)) (42) 
^^^ = (■^^^^•■^^(^0)) (43) 



Although Eqs. (41) and (42 1 both measure the correlations of the parallel velocity components, the former does so 
within the fixed lab frame, while the latter uses the moving body frame. Therefore, the first will be sensitive to changes 
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in the magnitude and direction of the velocity vector, while the second will only register changes in the magnitude. At 
short times, we can expect the decay in correlations of the parallel velocity components to be determined by the local 
configuration of the suspension (where the flow generated by the nearby particles can act to enhance or suppress the 
swimming motion), and at long times by the particle collisions, which will reorient the particle's swimming direction. 
This means the initial decay should be the same for both Cy and Cy, since the particle has not had time to reorient 
itself; at long times however, correlations measured within the lab frame of reference should decay to zero, whereas 
those measured within the body reference should reach a finite value (determined by the average swimming speed 
of the squirmers). For the decay in correlations of the perpendicular velocity components Cy we expect a behavior 
analogous to the short-time decay of the parallel components, since it is due entirely to the flow field generated by 
the neighboring squirmers. 

III. VALIDATION 
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FIG. 5. (color online) Swimming speed of an isolated squirmer in a periodically replicated cubic simulation box of length 
L/A = 128 at Re = 0.01 as a function of time (in simulation units), (a) Neutral swimmer at various particle sizes a/ A. (b) 
puller of size a/A = 6 for various swimming modes < a < 5. Velocities are scaled by the theoretical value for the swimming 
speed U = 2/3Bi. 

The first obvious test of our simulation method is to make sure that an isolated swimmer will move at the expected 
velocity U — 2/3Bi, regardless of the particle size or the value of the second squirming mode i?2- We performed 
simulations for a single squirmer, inside a periodically replicated cubic simulation box of dimension L = 128A (A 
is the grid spacing), for various particles sizes and squirming modes. Figure [s] shows the results obtained for a 
neutral squirmer (a — 0), for particle sizes a/A = 4,5,6,8, and for various pullers (a > 0), with a particle size of 
a/ A = 6. The particle Reynolds number Re — pUa/rj and the width of the diffuse interface C/A were the same for 
all the simulations, 0.01 and 2, respectively. In all cases, the particle's (average) swimming velocity shows excellent 
agreement with the theoretical predictions. For the neutral squirmers, the swimming velocity is within ~ 2% of the 
exact value, regardless of the particle diameter, although the speed shows a small decrease with increasing particle 
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size. The small oscillations exhibited by the velocity are due to discretization errors, as the number of grid points on 
the two hemispheres of the particle surface will vary depending on the relative position of the particle center within 
the computational bins. Similar agreement is obtained for the pullers (pushers), although the discretization error 
increases with increasing \a\ (at fixed particle size a/A). As such, for the system sizes we have considered, we are 
limited to moderate values of |a|< 5. 
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FIG. 6. (color online) Azymuthaly averaged steady-state fluid velocities for a single puller a = +2 of size a/A = 6 swimming 
along the £-axis, within a periodically replicated cubic simulation box of size L/A = 128. (a) Fluid velocity vectors within the 
laboratory frame, the red and blue arrows show the simulation and analytical results, respectively, (b) Fluid velocity stream 
lines within the laboratory frame, the colored lines (color-coded with respect to the magnitude of the velocity) represent the 
analytical solution, while the gray lines show the simulation results, (c) Same as Figure (b) but within a reference frame moving 
with the particle. Due to discretization errors, streamlines can begin/end at the fluid-particle interface. Length scales have 
been scaled by the particle radius. 



To verify that the velocity field generated by our SP squirmers is correct, we compare with the exact (at Re = 0) 
analytical expression given in Eq. (|3|. The steady-state velocity fields generated by a single squirmer (L/A = 128, 
v4/A = 6, C/A = 2, a = +2, Rc = 0.01), along with the corresponding stream line plot (in the lab and particle 
reference frames), are shown in Figure pj Excellent agreement with the analytical results is obtained, although 
differences in the stream lines arise at large distances r/a ~ 4 for Q ~ ±7r/4 (with respect to the swimming axis 
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+z). Even though the fluid velocity within these regions is vanishingly small (compared to the swimming speed of 
the squirmer), a clear systematic deviation is observed in the direction of the stream lines. This is due to our use 
of periodic boundary conditions, which causes the stream lines to close in on themselves, in order to match at the 
boundaries of the simulation cell. A similar deviation is observed for the velocities along the swimming direction, but 
these occur at much larger distances. 
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FIG. 7. (color online) Perpendicular force felt by two (fixed) parallel squirmers {a = +2), of diameter a/ A — 10 and interface 
thickness ("/A = 2, as a function of distance r, for Re — 0.01 and a system size of L/A — 128. The dashed yellow line gives 
the fit to the expected functional form F/Fs — ^ln(r-/cr — 1) of the near-field force, with A = —0.327, while the dotted violet 
line gives the (exact) far-field force. All forces are given in units of the Stokes force for an inert sphere with the same speed 
{U = 2/3Bi). 



As in Ref. |351 we consider the interactions between two fixed squirmers at a distance r from each other, with 
parallel orientations. Figure [7] shows the results we have obtained for the force parallel to the displacement vector 
between two pullers (perpendicular to their swimming axes) , normalized by the Stokes force Fg — QurjaU for an inert 
particle moving with the same velocity [U = 2/3Bi). The simulations were carried out at Re = 0.01 for a box size 
of L/A = 128, with a particle radius of a/A = 6, and a swimming mode a = +2. The functional form for this 
perpendicular force has been given by Ishikawa et al.[55] 

Fnear CX log (e) (44) 

where e = r — a is the minimum separation distance between the surface of the particles (with a — 2a the particle 
diameter). Additionally, they also obtained exact far-field expressions for the force between the two squirmers from 
a generalization of Faxen's laws. 

At short distances the repulsive force experienced by the two squirmers is seen to follow the expected scaling rela- 
tionship up to r < 1.5cr, while the far-field force is approached asymptotically for r > Sa. The force is observed to be 
proportional to the swimming mode, such that F on a and F{a) = —F{—a). 



IV. RESULTS AND DISCUSSIONS 



A. Diffusion Coefficients 



We have studied the diffusive behavior of a semi-dilute suspension of identical non-buoyant squirmers, for various 
concentrations </> < 0.124 (with (f> — Aira^N/SV the packing fraction) and squirming modes a = 0, ±1, ±2. We work in 
reduced units, in which the density and viscosity of the host fiuid are unity p = /i = 1. All simulations were performed 
for squirmers of radius a/ A = 5 (interface width (/A = 2), in a cubic simulation box of length L/A = 64, for a 
particle Reynolds number of Re = pUa/fi ~ 0.05 {U is the swimming speed of an isolated squirmer). Although our 
systems are not very large L/a ~ 12, finite size effects have been shown to be small[2Bj, and since we will be mainly 
focused on studying the short-range hydrodynamic interactions among particles, they can be safely ignored. In what 
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follows, all quantities are presented in non-dimensionalized form, using as characteristic length, speed, and time units 
the particle radius a, swimming speed U = 2/3i3i, and the time required for an isolated squirmer to move a distance 
equal to its radius T = a/U. 
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FIG. 8. (color online) Translational and rotational diffusivities for a system of pullers (a = +2) at various concentrations, 
(a) Translational diffusivities obtained from the mean-squared displacements in the fixed-lab reference frame Dxit), and the 
body frame Drit) (with the displacement due to the intrinsic swimming motion suitably removed), dashed and solid lines, 
respectively; and (b) the rotational diffusivities DR{t). Lines of slope one, corresponding to purely ballistic motion, have been 
drawn for comparison. The arrows show the increase/decrease of the various diffusivities as a function of concentration <j). 



The diffusivities Dxit) and Dxit), defined in Eqs. (29 1 and (351, for a system of pullers (a = -1-2), at various 
concentrations, are shown in Figure |8] The standard diffusivities Dxit) show a ballistic regime which extends to 
very long times t ~ 100, after which the slope starts to decrease, indicating a transition towards a diffusive regime. 
However, purely diffusive motion (slope of zero) is only obtained for the highest concentrations and alpha values. This 
behavior has been analyzed in detail in ref. [28l In contrast, Drit) reaches the diffusive regime at much shorter times, 
for all concentrations and all non-zero values of alpha. We note that both quantities are measuring related phenomena, 
but the latter does so directly, since the motion due to the inherent particle swimming has been removed from the 
analysis. As a function of concentration, Dt shows a clear increase, which is consistent with the interpretation of 
the diffusive motion being caused by interactions with neighboring particles: as the concentration increases, so does 
the number of neighbors, and thus the strength of the interactions. In contrast, we see that Dt decreases with 
concentration. The reason for this is simple, since mainly measures the effect of the particles' own swimming, an 
increase in the diffusive behavior can only hinder this motion. The rotational diffusivities Db,, also shown in Figure |8) 
present the same basic features as D^, the onset of the diffusive regime is obtained within the same time interval 
and they exhibit a similar concentration dependence. Apart from a difference in scale, there is no clear differentiating 
factor between these two quantities. Similar results are obtained for all other non-zero values of a considered; where, 
as was pointed out by Ishikawa and Pedley[28], the effect of increasing (decreasing) a is the same as that of increasing 
(decreasing) the concentration, at least with respect to the diffusive motion of the squirmers, since they both lead to 
an increase (decrease) in the strength of the hydrodynamic interactions among particles. 
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B. Velocity fluctuations and correlations 




FIG. 9. (color online) Single-particle velocity time-correlation functions for a system of pul lers (a = -1-2) at (p = 0.124. The two- 
time relaxation approximation to the velocity correlation function c{^' (t)/Cv (0) (Eq. 461, with the amplitudes and relaxation 
times obtained from the simulation, is also shown (dash-dot red line). 

To have a better understanding of the how the diffusive motion arises, we now consider the single particle velocity 



correlation functions defined in Eqs.(39) - (|43[ ). We present in detail the results obtained for a system of pullers 
(a = -1-2) at (j) = 0.124, shown in Figure l9j but the same analysis applies for the other values of cj) (a) we have 
studied. The standard definition of the velocity auto-correlation Cy (t) results in a function that shows a two-time 
decay process, as evidenced by the shoulder that appears around 1 < t < 10. The initial decay (t ~ 1) is due to 
the changing configuration of neighboring particles (and their self-generated flow field), while the second, long time 
decay {t ~ 10) arises due to the reorientation of the particles' swimming direction, which is a consequence of the 
particle-particle collisions. Confirmation for this interpretation is provided by the three other correlation functions 
we have defined. Consider first the decay in correlations of the parallel velocity components, Cy(t) (lab frame) and 

Cy (particle frame). Both functions show the same initial decay, caused by local fluctuations in the velocity field, as 
the particle has not experienced any significant change in its swimming direction. However, at longer times, these 
two functions show a completely different behavior. The parallel correlations measured within the fixed lab system 
Cy{t) exhibit a slow (long-time) decay before eventually collapsing onto the full correlation function Cvit); while 
the correlations measured within the particle's frame of reference reach a plateau after the initial decay. This means 
that while the particle is constantly reacting to the flow field generated by its neighbors, by changing it's swimming 
speed, these fluctuations decay very fast, typically within the time it takes for the particle to travel its diameter (the 
characteristic length over which the flow field varies significantly). As such, the long-time decay in both Cv{t) and 
Cy{t) is due primarily to changes in the orientation of the swimming direction, and not to changes in the magnitude 
of the swimming speed. As expected, the correlations in the perpendicular velocities Cy{t) show a very fast decay, 
over a characteristic time equal to that of the initial short-time decay of the parallel velocity components. We also see 
a clear correspondence between the decay in correlations of the angular velocity Cq (t) and the perpendicular velocity 

Finally, the characteristic times and r/ for these two processes, governed by the hydrodynamic interactions and 
the particle collisions, can be obtained by assuming an exponential decay exp(— i/r) for the appropriate correlation 

function: Cy{t) and Cy{t) for and n, respectively (see Figure [9|. The two-time scale relaxation of the velocity 
fluctuations can then be approximated as 

d^\t;^,a) = x.(0,a)e-*/^-(*'") +Xi(0,a)e-*/^'(^'") (46) 
Xs = {SV ■ SV) 

XI = ul 

with C/|| = {V ■ 63) and where we have explicitly shown the dependence of the amplitudes x decay times t on 
the concentration and swimming mode of the particles. Figure [9] shows excellent agreement between this two-time 
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relaxation process (with all four parameters obtained from the simulation data) and the full correlation function Cy(i) 



defined in Eq (39). Furthermore, in what follows we show how the long-time decay can be simplified and interpreted 



in terms of a binary collision process between the swimmers. 

C. Scaling 

To confirm the interpretation we have given for the diffusive motion of the squirmers, and the emergence of the two 
time-scales, we perform a scaling analysis for the relevant parameters (diffusion coefficients and collision time scales), 
as a function of concentration (p. Although a similar analysis has been presented in ref. by Ishikawa and Pedley, 
they have not considered the effective diffusion Dt, which means that only the long-time behavior of the system, given 
by the particle collisions, was analyzed (the hydrodynamic interactions which give rise to the rapidly decaying velocity 
fluctuations are completely masked by the swimming motion). In a subsequent study '(JI], the authors considered the 
diffusion of fluid particles, as well as inert colloidal particles, in a suspension of swimmers, and in this case they 
were able to observe the emergence of a second (shorter) time-scale. Our results, which are in agreement with their 
scaling arguments, provide a complementary view of the diffusive motion of these squirmer systems. The benefit of 
the analysis we propose lies in the fact that all relevant time and length scales can be obtained just from the motion of 
the squirmers, i.e., there is no need to consider the motion of the fluid or to introduce inert (non-swimming) particles. 
This results in simpler simulations and will also be relevant when trying to compare with experimental data. 

1. Correlation times 

First, lest us consider the long-time dynamics of the system, which is determined by the particle-particle collisions 
of the swimmers. From the kinetic theory of gases[53', we know that the mean free path A (or average distance 
between collisions) should be proportional to the inverse of the concentration of particles. In reduced units, we have 

A=^(^.^<^)~' (47) 

with (j) the packing fraction, and CTc the "collision" or cross-section diameter, which is not necessarily equal to the 
actual physical diameter of the particle. This difference is due to the self-generated flow profile around each particle, 
and as such should only depend on the swimming mode a, not on the concentration of particles. For the moment 
though, let us assume that we are dealing with hard-spheres Uc = 2, the average time between collisions is simply 
t^^ X^^/Uc, which gives 

= ^(C/.<^)-^ (48) 

where Uc is the average velocity of the particles. In thermal systems, this velocity would be determined by the 
reservoir and be independent of concentration, so that one obtains the well known (/)~^ dependence for the collision 
time t^^ oc (p~^ . Although the velocity of our squirmers is not concentration independent, to first order, we can safely 
assume that the particles are all swimming at an average velocity which is near the swimming velocity of an isolated 
squirmer, Uc — 1. Therefore, the collision time for our systems should show the same concentration dependence 
predicted by kinetic theory, tc oc (f)~^. To obtain the collision times tc from the decay times r computed from the 
simulations, we use the Enskog approximation [5^ [55] 

tc = 2t/3 (49) 



Our results, shown in Figure 10 confirm the scaling predictions for the long-time decay t[. = 2r;/3. Of special interest 
is the difference between for the different squirming modes, with the collision time decreasing with increasing a. This 
can be explained by an increase in the effective size of the particle, as the mean free path A is inversely proportional 
to the collision diameter ac- We also show the mean collision time t^^ for an equivalent system of hard-spheres. The 
fact that the collision times for the squirmers are significantly larger than the corresponding hard-spheres values is 
surprising, since it implies that the squirmers have a cross-section diameter which is smaller than the hard-sphere 
diameter of the particle [tc cx cr~^). The reason for this lies in the collision dynamics of the squirmers. 

The scaling of the fast time-scale (or t*) is harder to elucidate, since it shows only a very weak concentration 
dependence. Given that is related to the velocity fluctuations caused by the short-range hydrodynamic interactions 
between particles, it should be related to the time it takes for the flow region around a particle to change: this is 
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FIG. 10. (color online) (a) Scaling of the "collision" times, t^. and t", with concentration (f)- Also shown is the collision time t^^ 
for an equivalent system of hard spheres, with a (scaled) velocity of Uc = 1, given by the kinetic theory of gases, (b) Scaling of 
the translational Dt and rotational Dr diffusion coefficients with concentration, (c) Scaling of the velocity fluctuations with 
concentration. 



essentially the time necessary for a particle to swim its diameter t ~ 2. Our results, also shown in Figure [TOj agree 
with this rough estimate and also indicate a power-law behavior with an exponent of ~ —1/6. Although the decrease 
with concentration seems clear, it is very difficult to accurately measure such small variations and we have no suitable 
explanation for the value of this exponent (given by a fit to the data). These results are analogous to those obtained 
by Ishikawa et al. |31| for the time-scale of diffusing fluid particles in a suspension of squirmers. 
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2. Diffusion coefficients 

For the scaling of the translational diffusion coefficients, a simple dimensional analysis yields 

Dt oc C/2 tc (50) 
where Uc and tc are the characteristic velocity and time scales of the "collision" process which gives rise to the diffusive 



mot ion [31]. For the standard diffusion coefRcient Dt (Eq. (36 1), the reorientation of the particles is caused by the 
particle-particle collisions, thus tc = t\ and Dt oc (f)~^ (assuming a constant swimming velocity). This behavior has 
been analyzed in detail in ref. [28l More relevant to our study is the scaling of the effective diffusion coefficient Dt 



(Eq. (38 1), in which the motion due to the particle swimming has been removed. Our results, shown in Figure 10 
indicate a linear dependence with concentration Dt on (p. This is precisely the scaling behavior reported by Ishikawa 
et al.|31j for the diffusion of inert/fluid particles, and the same concentration dependence that has been observed 
experimentally for the diffusion of tracer particles in swimming suspensions (for both pushers and pullers) \29\ 130) . 
This is not surprising, since we can consider that this diffusive motion arises due to collisions between the particle 
(swimmer or not) and the fluid. By definition, it is clear that the velocity and time scales cannot be the same 
swimming speed or collision time used to define Dt- In this case, the characteristic velocity scale will be set by the 
velocity fiuctuations Sv around the average swimming velocity, which can be directly measured by computing the 
velocity components perpendicular to the swimming axis of the particles V-^. For these fluctuations, our simulations 



indicate a square-root dependence with concentration Sv cx (see Figure 10). The characteristic time-scale for this 
process will be the time necessary for the fluid flow surrounding a given particle to show considerable variations, and 
this can be estimated by the time it takes the particle to travel its own diameter tc cx 1. We thus recover the linear 
dependence given by the simulations. The same scaling behavior is obtained for the rotational diffusion coefficient 
Dr (X (f>. 

3. Collision diameters 

So far, we have not considered the effect of the squirming mode a on the diffusive properties of the swimmers. While 
it is clear that increasing the magnitude of a will give rise to stronger hydrodynamic interactions, and thus increase 



the diffusion of the particles (as can be seen in Figure 10), it is not clear how this dependence can be quantified. We 
propose that a, which defines the strength/range of the self-generated fiuid-flow around a squirmer (Eq. ([3])), can be 
directly related to the effective "collision" diameter of the swimmers. From the expression of the mean free path. 



Eq. (47), we obtain the following for the collision radius Vc — (Tc/2 



= (3V2UctcA (51) 



Assuming tc = /(a) <l> ^, with f{a) a function only of the squirming mode a, and considering the swimming velocity 
to be constant {Uc = 1), we arrive at a concentration independent collision radius 

rc'^[f{a)r'/^ (52) 
As f{a) is a decreasing function of a, the effective collision radius of the particle should increase with the magnitude 



of the squirming mode. Our results for the colfision radius, given by Eq. (51 ), with Uc = C/y = (V • 63) and tc = 2/3t' 
obtained directly from the simulations, are shown in Figure |11| Allowing for the large uncertainties in measuring 
the decay times at low concentrations, we obtain very good agreement with the previous scaling analysis: the radii 
show only a small variation with 0, and a clear distinction is observed as \a\ is increased. We note however that the 
pushers seem to present a larger collision radius than the pullers, something which is difficult to explain with the 
simple collision model we have presented. 

Also shown in Figure[Tl]is the mean collision radius Vc as a function of q. Although more data is needed to accurately 
specify the functional form of /(a), a simple linear dependence predicts a value of rda = 0) which is an order of 
magnitude smaller than the hard-sphere radius of the particles rc — 0.1. This is consistent with our simulations, 
from which we were unable to obtain reliable estimates for r/, because the velocity correlation functions exhibited 
no substantial decay over the time scales we studied t ~ 10^ ~ 10^. However, we were able to obtain estimates for 
the effective diffusion coefficients Dt of these neutral squirmers, which are an order of magnitude smaller than the 
corresponding values for a — ±1. This means that the difference in collision times, with respect to a corresponding 



system of hard-spheres (see Figure 10), will be even larger for the neutral squirmers than it is for the swimmers with 



non-zero a values. This is consistent with an interpretation in terms of hard-spheres when analyzing the differences 
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FIG. 11. (color online) Effective collision radius of the squirmers as a function of concentration Results were obtained by 
using the mean free path expression for a system of hard-spheres given by the kinetic theory of gases, with the collision time 
given by t\ — r;/16 (r; is the decay time in the velocity correlations given by particle-particle collisions). 



among squirmers (larger a equivalent to larger collision radius Vc), but it can appear contradictory when comparing 
to actual hard-sphere systems. After all, we could expect the systems for a = 0, in which the flow profile due to the 
squirming motion is most localized {U{r) cx r~^), to be closest to a system of hard-spheres. Yet our results indicate 
that the collision radius of the neutral squirmers is an order of magnitude smaller than the actual radius of the 
particles. As we have already mentioned, this is due to the difii'erence in the collision dynamics of the particles, which 
differ dramatically from that of hard-spheres [33l l54]. Although the exact collision process will depend not only on 
the relative velocity of the particles, but also on their relative orientations, the squirming motion results in deflection 
angles which are (on average) smaller than the corresponding hard-sphere values. In essence, this means that the flow 
fields generated by the particles allows them to swim past each other with a relatively small change in their swimming 
direction (again, as compared to with an actual hard-sphere collision). 



Finally, using Eq. (51 1, the velocity correlation function for a suspension of squirmers (Eq. (46)) can be reduced to 
the following functional form 



(53) 



where the long-time decay process is expressed only in terms of the average swimming speed C/|j((/), a) and the 
hard-sphere collision diameter ada) of the squirmers. The former has a weak linear dependence on concentration 
?7|| = 1 — A{a)(j), while the latter is concentration independent, as was shown in Fig. 11 The scaling behavior 
for the short time process is still not understood, but our simulation results suggest a power law behavior for the 
concentration dependence of the amplitude Xs and decay time of the velocity fluctuations, with an exponent of 1/2 
and —1/6, respectively; i.e., Xs cx and Ts oc A more detailed analysis is required to firmly establish these 

scaling relationships, as well as to determine the dependence on the swimming parameter a. Work along these lines 
is in progress. 
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V. CONCLUSIONS 



We have investigated the hydrodynamic interactions of suspensions of squirmers using a modified version of the 
smoothed profile method (SP) for particle dispersions. The SP method allows one to fully resolve the hydrodynamic 
interactions in many particle dispersions in an accurate and efficient manner, and we have shown how it can be 
extended to systems with slip boundary conditions, such that it is possible to describe squirmers (active swimmers 
which move due to self-generated surface tangential velocities) . The validity of the method was confirmed by comparing 
the simulation data with exact results for the case of a single swimmer, for which we recover the correct swimming 
speed and are able to accurately reproduce the fluid flow generated by the squirming motion, and for two aligned 
swimmers at a fixed distance, for which we recover the expected hydrodynamic force. The advantage of the SP 
method for swimming particles, as opposed to Stokesian Dynamics (which has been successfully, and extensively, used 
to study these systems) is its applicability to particle dispersions in complex fluids. This is relevant in the case of 
swimming micro-organisms, since the role of the nutrient and the possibility of having a non-Newtonian host fluid 
must be considered when comparing with experiments. 

In this paper we have analyzed the effect of the hydrodynamic interactions on the motion of semi-dilute squirmer 
suspensions, up to volume fractions of (/) = 0.124, for various swimming modes |a| < 2. Although we have no yet 
included thermal fluctuations in our description, the swimming motion of the particles gives rise, over sufficiently 
long time scales, to a diffusive regime. In order to distinguish between the contributions due to the hydrodynamic 
interactions, caused by the squirming motion, and those due to the particle-particle collisions, which are the two 
basic mechanisms responsible for the diffusive motion, we have analyzed the particle dynamics in terms of movement 
due to the inherent swimming of the particles, and that due to the (hydrodynamic) interactions among them. This 
is easily done by looking at the motion from the particle's own frame of reference, i.e., decomposing the motion 
parallel and perpendicular to the swimming axis. This analysis has allowed us to demonstrate the appearance of two 
distinct time scales within our system, one related to the time between particle-particle collisions, the other to the 
fluid-particle interactions. This two-time scale nature of the particle interactions in swimming suspensions can be 
clearly seen in the two-time relaxation of the velocity correlation function. We are thus able to define an effective 
hydrodynamic diffusion coefficient (corresponding to the short-time fluid/particle interactions), in which the self- 
motion of the particle has been removed, which shows a linear scaling with concentration. In contrast, the standard 
diffusion coefficient is inversely proportional to the concentration of swimmers. This is in agreement with simulation 
and experimental results on tracer diffusion in swimming suspensions. Additionally, since the long-time dynamics 
of the system is related to the particle-particle collisions, we have used the well-known results from kinetic theory 
to deduce an effective, concentration independent, collision radius for our swimmers. Due to the complex collision 
dynamics of these particles, this collision radius is actually smaller than the hard-sphere radius of the particle, and 
increases with increasing a. 
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